%\chapter{FDTD Modelling of Metamaterials}
\chapter{Introduction to Finite Difference Time--Domain Technique}
\section{The Yee Algorithm}
The original algorithm was\index{Yee algorithm} proposed by K. S. Yee\index{Yee, K. S.} in 1966~\cite{Yee1966}. The derivatives in Faraday's Law\index{Faraday's Law} and Ampere's Law\index{Ampere's Law} are replaced with difference equations. The whole problem space is discretised and divided into cells such that electric and magnetic fields are staggered at half spatial steps in unit cell. By advancing the simulation in time, future values of electric and magnetic field components are computed using past values. This chapter gives provides a brief treatment of FDTD\index{FDTD} in 1D and 2D. Wave propagation in free space and with a homogeneous scatterer are discussed using conventional and dispersive FDTD\index{FDTD!dispersive} methods.
\section{FDTD Update Equations in One Dimension (1D)}
Faraday's Law\index{Faraday's Law} and Ampere's Law\index{Ampere's Law} in differential form are given by
\begin{equation}
\centering
\nabla \times \textbf{E} = - \dfrac{\partial \textbf{B}}{\partial t},
\label{Faraday's Law}
\end{equation}
\begin{equation}
\centering
\nabla \times \textbf{H} = \dfrac{\partial \textbf{D}}{\partial t}.
\label{Ampere's Law}
\end{equation}
Assume electric field only has $x$--component, that is $E_x$, and magnetic field only has $y$--component $H_y$. The curl of electric field can be expanded as
\begin{equation}
\centering
\nabla \times \textbf{E} = - \mu \dfrac{\partial H_y}{\partial t} \hat{y} = \left| \begin{array}{ccc} \hat{x} & \hat{y} & \hat{z} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ E_x & 0 & 0 \end{array} \right| = \dfrac{\partial E_x}{\partial z} \hat{y}.
\label{Faraday's Law-Curl-Expansion}
\end{equation}
Similarly, the curl of magnetic field can be written as
\begin{equation}
\centering
\nabla \times \textbf{H} = \epsilon \dfrac{\partial E_x}{\partial t} \hat{x} = \left| \begin{array}{ccc} \hat{x} & \hat{y} & \hat{z} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ 0 & H_y & 0 \end{array} \right| = - \dfrac{\partial H_y}{\partial z} \hat{x}.
\label{Ampere's Law-Curl-Expansion}
\end{equation}
Rearranging terms,
\begin{equation}
\centering
\dfrac{\partial H_y}{\partial t} = - \dfrac{1}{\mu} \dfrac{\partial E_x}{\partial z},
\label{Faraday's Law-Curl-Expansion-1D}
\end{equation}
\begin{equation}
\centering
\dfrac{\partial E_x}{\partial t} = - \dfrac{1}{\epsilon} \dfrac{\partial H_y}{\partial z}.
\label{Ampere's Law-Curl-Expansion-1D}
\end{equation}
These two scalar equations drive the FDTD\index{FDTD} algorithm. First, magnetic field is calculated using equation~\ref{Faraday's Law-Curl-Expansion-1D} and then electric field (equation~\ref{Ampere's Law-Curl-Expansion-1D}) is calculated from the magnetic field obtained in~\ref{Faraday's Law-Curl-Expansion-1D}. Electric and magnetic fields are staggered at half time steps. Similarly, magnetic and electric fields in space are also spaced at half spatial steps. The algorithm is normally referred to as a leap--frog scheme\index{leap--frog}.
\begin{figure}[here]
\centering
\begin{tikzpicture}[xscale=1.2,yscale=1.2,font=\small]
	\def\XD{0cm}
	\def\YD{0cm}
	\coordinate[label=left:$E^{n=0}$] (E0) at (\XD,\YD+0.3cm);
	\foreach \x in {0.0cm, 1.0cm, 2.0cm, 3.0cm, 4.0cm, 5.0cm, 6.0cm, 7.0cm}
		\draw[thick] (\XD+\x, \YD+0cm) rectangle (\XD+\x+1.0cm, \YD+0.6cm);
	\foreach \x/\y in {0.0cm/0, 1.0cm/1, 2.0cm/2, 3.0cm/3, 4.0cm/4, 5.0cm/5, 6.0cm/6, 7.0cm/7}
		\draw (\XD+\x+0.5cm,\YD+0.3cm) node {$E[\y]$};
	
	\def\XD{0.5cm}
	\def\YD{1.5cm}
	\coordinate[label=left:$H^{n=\frac{1}{2}}$] (H1/2) at (\XD,\YD+0.3cm);
	\foreach \x in {0.0cm, 1.0cm, 2.0cm, 3.0cm, 4.0cm, 5.0cm, 6.0cm, 7.0cm}
		\draw[thick] (\XD+\x, \YD+0cm) rectangle (\XD+\x+1.0cm, \YD+0.6cm);
	\foreach \x/\y in {0.0cm/\frac{1}{2}, 1.0cm/1\frac{1}{2}, 2.0cm/2\frac{1}{2}, 3.0cm/3\frac{1}{2}, 4.0cm/4\frac{1}{2}, 5.0cm/5\frac{1}{2}, 6.0cm/6\frac{1}{2}, 7.0cm/7\frac{1}{2}}
		\draw (\XD+\x+0.5cm,\YD+0.3cm) node {$H[\y]$};

	\draw[dashed] (-2.5cm,2.5cm) -- (9.0cm,2.5cm);

	\def\XD{0cm}
	\def\YD{3cm}
	\coordinate[label=left:$E^{n=1}$] (E1) at (\XD,\YD+0.3cm);
	\foreach \x in {0.0cm, 1.0cm, 2.0cm, 3.0cm, 4.0cm, 5.0cm, 6.0cm, 7.0cm}
		\draw[thick] (\XD+\x, \YD+0cm) rectangle (\XD+\x+1.0cm, \YD+0.6cm);
	\foreach \x/\y in {0.0cm/0, 1.0cm/1, 2.0cm/2, 3.0cm/3, 4.0cm/4, 5.0cm/5, 6.0cm/6, 7.0cm/7}
		\draw (\XD+\x+0.5cm,\YD+0.3cm) node {$E[\y]$};

	\def\XD{0.5cm}
	\def\YD{4.5cm}
	\coordinate[label=left:$H^{n=1\frac{1}{2}}$] (H11/2) at (\XD,\YD+0.3cm);
	\foreach \x in {0.0cm, 1.0cm, 2.0cm, 3.0cm, 4.0cm, 5.0cm, 6.0cm, 7.0cm}
		\draw[thick] (\XD+\x, \YD+0cm) rectangle (\XD+\x+1.0cm, \YD+0.6cm);
	\foreach \x/\y in {0.0cm/\frac{1}{2}, 1.0cm/1\frac{1}{2}, 2.0cm/2\frac{1}{2}, 3.0cm/3\frac{1}{2}, 4.0cm/4\frac{1}{2}, 5.0cm/5\frac{1}{2}, 6.0cm/6\frac{1}{2}, 7.0cm/7\frac{1}{2}}
		\draw (\XD+\x+0.5cm,\YD+0.3cm) node {$H[\y]$};


	\draw[thick, ->, >=stealth] (-1.5cm, -0.5cm) -- (-1.5cm, 5.5cm);
	\coordinate[label=above:$time$] (time) at (-1.5cm,5.5cm);
	\coordinate[label=left:\textbf{Past}] (Past) at (-1.5cm,2.2cm);
	\coordinate[label=left:\textbf{Future}] (Future) at (-1.5cm,2.8cm);

\end{tikzpicture}
\caption{FDTD leap--frog\index{leap--frog} scheme to update future fields from past fields}
\label{FDTD-LeapFrog}
\end{figure}
By examining equations~\ref{Faraday's Law-Curl-Expansion-1D} and \ref{Ampere's Law-Curl-Expansion-1D} it can be intuitively observed that a wave with $E_x$ and $H_y$ components will propagate in $\hat{z}$ direction and space derivatives with respect to $z$ indicate variation of electric and magnetic fields in space along $z$--axis.

These equations need to be discretised before they can be implemented as a computer simulation. Both the time and space derivatives are discretised to obtain difference equations. The difference equations for space and time derivatives use central difference\index{central difference} scheme and expressed as
\begin{equation}
\centering
\dfrac{\partial F^n(x)}{\partial t} \equiv \dfrac{F^{n+\frac{1}{2}}(x)-F^{n-\frac{1}{2}}(x)}{\Delta t},
\label{central-difference-time}
\end{equation}
\begin{equation}
\centering
\dfrac{\partial F^n(x)}{\partial x} \equiv \dfrac{F^n(x+\frac{1}{2})-F^n(x-\frac{1}{2})}{\Delta x}.
\label{central-difference-space}
\end{equation}
Following the treatment in chapter 3 of~\cite{JBSchneiderUFDTD}, electric and magnetic fields as functions of time and space can be written as
\begin{equation}
\centering
H_y (z,t)=H_y (k\Delta _z,n\Delta _t)=H^n_y[k],
\label{Hy-of-z_and_t}
\end{equation}
\begin{equation}
\centering
E_x (z,t)=E_x (k\Delta _z,n\Delta _t)=E^n_x[k].
\label{Ex-of-z_and_t}
\end{equation}
The final equations for electric and magnetic fields driving the FDTD algorithm are\index{FDTD!update equations}
\begin{equation}
\centering
H^{n+\frac{1}{2}}_y \left[k+\frac{1}{2}\right]=H^{n-\frac{1}{2}}_y \left[k+\frac{1}{2}\right] + \dfrac{\Delta t}{\mu \Delta x} \left( E^{n}_x \left[k\right] - E^{n}_x \left[k+1\right] \right),
\label{Hy-1D-Simple-FDTD-Driver}
\end{equation}
\begin{equation}
\centering
E^{n+1}_x \left[k\right]=E^{n}_x \left[k\right] + \dfrac{\Delta t}{\epsilon \Delta x} \left( H^{n+\frac{1}{2}}_y \left[k-\frac{1}{2}\right] - H^{n+\frac{1}{2}}_y \left[k+\frac{1}{2}\right] \right).
\label{Ex-1D-Simple-FDTD-Driver}
\end{equation}
The FDTD\index{FDTD} code for this 1D simulation is given in appendix \ref{App:1DFDTDFreeSpaceMatlab} and uses first order absorbing boundary condition (ABC)\index{boundary condition!absorbing (ABC)}.
\section{Stability and Accuracy}
\subsection{Courant Stability Criterion}
Stability of the FDTD\index{FDTD} method is governed by the temporal time step $\Delta t$. If it is too large then the simulation would become unstable. In FDTD\index{FDTD} the fields exist at discrete locations. In one iteration or temporal step, energy should not propagate any further than the neighbouring field or node. This is known as the Courant Stability Criterion\index{Courant stability criterion} and is expressed as
\begin{equation}
\centering
c \Delta t \leq \left[ \dfrac{1}{\Delta x^2} + \dfrac{1}{\Delta y^2} + \dfrac{1}{\Delta z^2} \right]^{-\frac{1}{2}}.
\label{Courant-Stability-Criterion}
\end{equation}
For the case where spatial differentials are equal (denoted by $\delta$), this reduces to
\begin{equation}
\centering
S_c = \dfrac{c \Delta t}{\delta} \leq 1.
\label{Courant-Stability-Criterion-Equal-Differentials}
\end{equation}
Here $S_c$ is known as the Courant number\index{Courant number ($S_c$)}. It should be less than unity in order for a stable solution.
\begin{figure}[htb]
\centering
\mbox{\subfigure[$S_c = 1$]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Stable1DFDTD.pdf}}\quad\subfigure[$S_c = 1.1$]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Unstable1DFDTD.pdf}}}
\caption{$S_c$ and stability of FDTD}
\label{Unstable-FDTD}
\end{figure}
\subsection{Courant Number as Scaling Factor}
In one time step ($\Delta t$) the amount of distance travelled by an electromagnetic wave (in cells) is equal to Courant number. Since $S_c$ cannot be greater than unity, an EM wave can only travel one cell or $\Delta$ during one time stepping of algorithm. Simulations with values of $S_c$ close to unity can rapidly become unstable in the presence of an obstacle or scatterer. Smaller values of $S_c$ will result in better stability but the trade--off is an increase in simulation time.
\subsection{Wavelength and Cell Size}
The size of wavelength should be considerably greater than minimum spatial step. Normally wavelength is taken to be at least ten times the size of a single cell~\cite{Taflove2000}. It makes sense that there shouldn't be too much fluctuations between two spatial steps and wavelength should be spread over a reasonable number of cells. Figure~\ref{Pulsewidth-comparison} shows different sizes of Gaussian pulses\index{source!Gaussian pulse} in a 1D simulation. Decreasing pulse width below 10 cells results in noticeable fluctuations and this can result in inaccurate results.
\begin{figure}[htb]
%\centerline{\includegraphics[]{GPGPU}}
\centering
\begin{tikzpicture}
	\draw[thick] (0cm, 0cm) rectangle (4cm, 4cm);
	\foreach \x in {5cm, 6cm, 7cm, 8cm}
		\draw[thick] (\x, 0cm) rectangle (\x + 1 cm, 4cm);
	\draw (0,3.5) cos (4,0.5);
	\draw (5,3.5) cos (9,0.5);
\end{tikzpicture}
\caption{Wavelength and cell size}
\label{Wavelength-Vs-Cell-Size}
\end{figure}
% Reference for double figure: http://www.johndcook.com/blog/2009/01/14/how-to-display-side-by-side-figurs-in-latex/
\begin{figure}[htb]
\centering
\mbox{\subfigure[Gaussian pulse of width 2]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Pulsewidth2.pdf}}\quad\subfigure[Gaussian pulse of width 4]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Pulsewidth4.pdf}}}\\
\mbox{\subfigure[Gaussian pulse of width 8]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Pulsewidth8.pdf}}\quad\subfigure[Gaussian pulse of width 16]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Pulsewidth16.pdf}}}
\caption{Comparison of pulse width}
\label{Pulsewidth-comparison}
\end{figure}
Inaccuracies may also result from inappropriate selection of spatial or time steps. Figure~\ref{Unstable-FDTD} shows that FDTD\index{FDTD} quickly becomes unstable if $S_c$ becomes greater than unity. Further to this, FDTD\index{FDTD} is inherently dispersive\index{FDTD!dispersion} in nature and a pulse will spread as it travels over a period of time. While these inaccuracies cannot be completely eliminated, they can be reduced with proper selection of simulation parameters.
\section{Grid Termination and Boundary Conditions}
This section is a discussion on a number of grid termination conditions also known as boundary conditions. A graphical picture of different boundary conditions is given in figure~\ref{Boundary-Conditions}.
\subsection{Perfect Electric Conductor (PEC)\index{boundary condition!perfect electric conductor (PEC)}}
Where electric field nodes are located at the outermost boundary these can be set to zero to simulate a PEC boundary. A wave incident on PEC will be completely reflected with a 180 degree phase shift.
\subsection{Perfect Magnetic Conductor (PMC)\index{boundary condition!perfect magnetic conductor (PMC)}}
PMC is modelled by setting magnetic field nodes at edges set to zero. A PMC boundary will reflect incoming wave without any change in phase.
\subsection{Periodic Boundary Condition (PBC)\index{boundary condition!periodic (PBC)}}
In PBC incident waves will wrap around the and appear on the opposite edge. It is used to model the domain as a circular or continuous space. PBC boundary is implemented by incorporating the fields on edges in the FDTD update equations\index{FDTD!update equations}.
\subsection{Absorbing Boundary Condition (ABC)\index{boundary condition!absorbing (ABC)}}
An ABC will completely absorb incident fields without any reflection. This is useful for antenna radiation and scattering problems. There are a number of techniques for modelling ABCs, most popular being perfectly matched layer or PML\index{boundary condition!PML}.
\begin{figure}[H]
\centering
\mbox{\subfigure[Incident wave]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Incident_Wave.pdf}}\quad\subfigure[Reflection from PMC]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_PMC_Reflected.pdf}}}
\mbox{\subfigure[Incident wave]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Incident_Wave.pdf}}\quad\subfigure[Reflection from PEC]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_PEC_Reflected.pdf}}}
\mbox{\subfigure[Incident wave]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Incident_Wave.pdf}}\quad\subfigure[Periodic boundary (PBC)]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_PBC.pdf}}}
\mbox{\subfigure[Incident wave]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_Incident_Wave.pdf}}\quad\subfigure[Absorbing boundary (ABC)]{\includegraphics[scale=0.375, trim=3.75cm 8.75cm 4.5cm 9.25cm, clip]{Figures/FigCh02_ABC.pdf}}}
\caption{Boundary Conditions: Incident wave is moving towards right.}
\label{Boundary-Conditions}
\end{figure}
